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Abstract 

Drug discovery today is a complex, expensive, and time-consuming process with high attrition rate. A more 
systematic approach is needed to combine innovative approaches in order to lead to more effective and efficient 
drug development. This article provides systematic mathematical analysis and dynamical modeling of drug effect 
under gene regulatory network contexts. A hybrid systems model, which merges together discrete and continuous 
dynamics into a single dynamical model, is proposed to study dynamics of the underlying regulatory network under 
drug perturbations. The major goal is to understand how the system changes when perturbed by drugs and give 
suggestions for better therapeutic interventions. A realistic periodic drug intake scenario is considered, drug 
pharmacokinetics and pharmacodynamics information being taken into account in the proposed hybrid systems 
model. Simulations are performed using MATLAB/SIMULINK to corroborate the analytical results. 
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Introduction 

The ultimate goal of drug therapy is to modulate the 
phenotypic behavior of cells by altering the behavior 
of the gene and protein components of the cell [1]. 
This approach is possible because the phenotypic behav- 
ior of the cell reflects the dynamics of the gene and 
protein-based regulatory network. When it comes to drug 
therapeutics and disease modeling, the major goal is to 
understand how the system changes when perturbed and 
how to modify the system to achieve a desired outcome. 
To understand and exploit the complicated mapping 
between genome and phenome, especially in the context 
of drug discovery, it is critical to evaluate the regulatory 
interactions between the genes and proteins that form the 
gene regulatory network (GRN). To date, the hope of the 
rapid translation of "genes to drugs" has foundered on the 
reality that disease biology is complex and drug develop- 
ment must be driven by insights into biological responses 
[2]. A systems approach is crucial for moving biology 
from a descriptive to a predictive science [3,4] . This calls 
for appropriate modeling to establish a functional under- 
standing of disease-drug interaction, in order to better 
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predict drug effects and make drug discovery a faster and 
more systematic process. 

Pharmacokinetics (PK) is the study of what the body 
does to the drug, i.e., the absorption, distribution, 
metabolism, and excretion of the drug, and pharmacody- 
namics (PD) seeks to study what the drug does to the body, 
A salient challenge is to link a drugs PK information with 
PD characteristics to provide a better understanding of 
the time course of drug effect (PK/PD) after drug admin- 
istration [5]. Modeling and simulation tools are required 
to integrate PK and PD data and optimize drug regimens. 

A salient problem is finding a dosing regimen of a 
drug candidate that is both efficacious and safe [6]. Tra- 
ditionally, drugs have been administered on an experi- 
mental basis, but it is virtually impossible to optimize 
dosing regimens using strictly empirical methods, espe- 
cially since different patients may respond differently to 
the same drug dosage [7] . Moreover, traditionally design- 
ing the dosing regimen to achieve some desired target 
goal such as relatively constant serum concentration may 
not be optimal because of underlying dynamic biological 
networks. For example, Shah et al. [8] demonstrate that 
BCR-ABL inhibitor dasatinib, which has greater potency 
and a short half-life, can achieve deep clinical remission 
in CML patients by achieving transient potent BCR-ABL 
inhibition, while traditionally approved tyrosine kinase 
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inhibitors usually have prolonged half-lives that result in 
continuous target inhibition. A similar study of whether 
short pulses of higher dose or persistent dosing with lower 
doses has the most favorable outcomes has been carried 
out by Amin et al [9] in the setup of inactivation of HER2- 
HER3 signaling. Finding an optimal dosing regimen based 
on the dynamics of biological systems and relevant PK/PD 
information is critically important. 

System modeling is emerging as a valuable tool in thera- 
peutics to address these challenges [3,10-12]. The process 
begins with building a quantitative model of a biological 
system. Consequences of particular perturbations, such as 
optimal dosing regimens, optimal drug targets, or combi- 
national therapy, can be simulated in time courses using 
such models. In this study, we propose a hybrid sys- 
tems model for GRNs and incorporate a drugs PK and 
PD information by using a state-space approach. We first 
study drug effect assuming the drug target to be a gene 
or protein in the proposed drug perturbation model using 
dynamical system theory, considering the case of peri- 
odic drug intake and analytically deriving the conditions 
for the drug to be effective. We extend the analysis to 
the 2-gene case and then to the case of a network with 
multiple coupled genes and positive feedback loops. Sim- 
ulations are performed using MATLAB/SIMULINK to 
supplement our analytical results. 

Model formulation 

While discrete modeling leaves out many details, con- 
tinuous modeling includes so many details that com- 
putational demands preclude their applications to many 
larger systems. Hybrid systems, which aim to merge ideas 
from both continuous and discrete modeling into one 
paradigm, are appealing for GRN modeling under drug 
perturbations because biological systems are naturally 
nonlinear, have highly varied regulatory requirements, 
and possess a wide range of control strategies for meeting 
their needs. While some simple, local, feedback control 
methods can provide sufficient regulation of many more- 
or-less continuous cellular processes, the regulation of 
discontinuous processes possessing the character of com- 
putational decision making requires more elaborate reg- 
ulatory methods [13]. In particular, some genes display 
regulation in a thresholded switch-like manner [14]. 

Hybrid systems include a broad space of models and sys- 
tems. Several hybrid systems models have been developed 
for biological networks. Some of these have been used 
to perform reachability analysis to elucidate biologically 
meaningful properties. For example, the Lac operon sys- 
tem has been well studied both experimentally and using 
continuous models [15,16]. A hybrid model and use of a 
reachability algorithm were validated by comparison with 
experimental data and continuous models [17]. Other bio- 
logical hybrid systems analyzed in similar ways include 



the Delta-Notch decision process [18,19], GRNs of car- 
bon starvation [20], and nutritional stress response [21] 
in Escherichia coli. As far as we know, the only hybrid 
systems modeling concerning treatment or drug effects is 
contained in our earlier work [22]. 

Gene regulation can be modeled by rate equations 
expressing the difference between rate of production and 
rate of degradation [23,24]. We adopt the general model 

*i =fi(x) - Yi*h (1) 

where X{ > 0 corresponds to the concentrations of 
proteins encoded by genes in the network and can be 
interpreted as the gene expression level. is a general 
nonlinear function and represents the rate of synthesis. 
It can be approximated by a sigmoidal function or a unit 
step function, and unit step function is used in this article. 
Yi%i is the rate of degradation. To use hybrid systems and 
incorporate drug effect, we propose the following model 
for a GRN of N genes under drug perturbation: 
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- YiXi + Pi^i ~ Yi &txu Vi = 1, 2, ... , N, 

where the last two terms on the right-hand side of 
Equation (2) correspond to drug perturbation u. /3 > 0 
and y > 0 are synthesis and degradation rates, respec- 
tively. Ki > 1 is an integer representing the number of 
activation/synthesis terms. Of. and E J - describe how other 
genes affect gene /. They are the functions of s + {xj 9 0^ +i ) 
and s~{Xj ) 0j~ i ) ) where is the unit step function, 

5~(.) = 1 — and the 0 terms are the corresponding 
threshold values. For each gene y, the set of threshold val- 
ues related to gene i is denoted by Tj, where £+* and t— i 

are indices of the threshold values, 0 < 6^ +i e Tj, and 

0 < Oj - ' e Tj. and O; represent the two sets of genes 
that affect the expression of gene i in different manners. 
Specifically, in this article, we consider 

^=5 + (^6>/ + 05-(^,6>/"0, (3) 

with E J . defined similarly. Q? t and E J . may be set to 0 or 
1, or different forms when appropriate threshold values 
are chosen. For example, Q\ = s + (#i,0~)s~(#i, oo) = 1 
and Qj = s+(x 2 , 0~)s~(x 2f 0^ ) = s-(x 2 ,0\). Qf and Ef 
describe how the drug u affect gene i. fif > 0 and y¥ > 0 
are the synthesis and degradation factors of the drug on 
gene i. PfQf and — yfEfxi are used when the drug is 
activating or repressing certain genes, respectively. Since 
most drugs are used to repress genes, only — y. M S^- is 
considered in the examples of this article. Note that y M is 
defined as a drug-effect factor, which is closely related to 
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the drug pharmacology model discussed in the following 
section. 

It should be kept in mind that the focus of this arti- 
cle is studying the effect of dosing, in particular, dosing 
regimens, on the expression of genes involved in a pathol- 
ogy by using hybrid systems theory. Whereas the sim- 
pler Equation (1) is widely accepted, it does not contain 
drug-effect terms. Equation (2) extends Equation (1) by 
including such terms. While the structure is intuitively 
reasonable and somewhat general, the actual details of the 
drug-effect terms are unknown. Finding the specific form 
of Equation (2) for a specific disease is a system identifi- 
cation problem, which is quite distinct from the analysis 
problem addressed in this article. We are addressing opti- 
mization of treatment intervention, given the system. The 
details of our analysis might change when the details of 
Equation (2) are clarified, but we expect that the hybrid 
systems approach taken in the article will go through with 
appropriate modifications in the mathematical details. 

We consider a 2-gene example to illustrate the feasibility 
of using hybrid systems for modeling drug effect. Specifi- 
cally, we assume that there are two interactive genes xi,X2 
that repress each other, and x\ is a disease gene which 
loses its self-regulation. We also assume that a drug targets 
X\ by reducing its expression level and providing a nega- 
tive feedback term —y^x\. The resulting 2-gene network 
is given by 



X\ 

x 2 



fos-ix^Ol) 



YiXi 

Y2*2 



(4) 
(5) 



where f$\, p2 are synthesis factors, 72 is a degradation fac- 
tor, and 0\, 0\ are threshold values, is a drug-effect fac- 
tor. Using dynamical systems theory, the state- trajectory 
schematic diagrams of this 2-gene network without and 
with drug input are obtained and plotted in Figures 1 and 
2, respectively. It is observed that without drug input, the 
gene expression level of x\ increases unbounded, while 
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Figure 1 State trajectory schematic of 2-gene example without 
drug intake. 




Figure 2 State trajectory schematic of 2-gene example with drug 
intake. 



with proper drug input, Pi/y" < Q\> the system converges 
to a new steady state, (fii/y^, P2/Y2)' 

We assume periodic drug intake and the drug concen- 
tration level in the effect-site follows exponential decay 
during each period X{, i.e., m(t) = ^ie~ kd ^~ kti \ where 
kxi < t < (k + 1)t/ and is the degradation factor. The 
response of gene expression levels of the two genes under 
periodic drug intake is shown in Figure 3. The state-space 
trajectory of gene expression level of x\ vs. the drug con- 
centration level u is given in Figure 4. A comparison of 
trajectory of the gene expression level x\ vs. X2 with and 
without drug are provided in Figures 5 and 6, respectively. 
It is observed that the drug is quite effective in bring- 
ing down the expression level of x\. The simulation study 
matches the theoretical analysis, as in Figures 1 and 2, that 
with proper drug intervention the system will converge to 
a new steady state, x\ — f$\/Y\ an d #2 = P2/Y2 — 1> while 
X2 — ► 0 and x\ —> 00 without drug input. 



drug concentration level 
x 1 (target gene) expression 

x 2 gene expression level 
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Figure 3 The state response under periodic drug intake. 
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Figure 4 The state-space trajectory under periodic drug intake. 

Parameter setting of Figures 3 and 4: X] (0) = 20, x 2 (0) = 0.7, r = 8, 
u(kr) = 24, ^ = 0.21, fa = fi 2 = 1, X2 = 1, 0, 2 = 10, 6>j = 2, 
A. d = 0.5. 



Pharmacology model 

The basis of clinical pharmacology is the fact that the 
intensities of many pharmacological effects are func- 
tions of the amount of drug in the body and, more 
specifically, the concentration of drug at the effect-site 
[5]. Historically, PK and PD were considered as sepa- 
rate disciplines; however, the information provided by 
these disciplines is limited if regarded in isolation [25]. A 
drug-effect factor y u is included in our proposed model 
(Equation 2), which is related to drugs PD characteristic 
(concentration-response) and its PK information (dose- 
concentration). In order to describe the time course of 
drug effect in response to different dosing regimens, the 
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Figure 5 The state-space trajectory with drug: t = 8, u(kr) = 24, 


q" = 0.21, Xd = 0.5. The rest parameter settings are the same with 


Figures 3 and 4. 
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Figure 6 The state-space trajectory without drug intake. The rest 


parameter settings are the same with Figures 3 and 4. 



integrated PK/PD model is indispensable because it builds 
the bridge between these two classical disciplines of phar- 
macology [25]. Following each dosing regimen, instead of 
a two-dimensional PK and PD relationship, the proposed 
approach enables a description of a three-dimensional 
dose-concentration-effect relationship. Specifically, PK 
and PD are linked through y u by a state-space approach to 
facilitate the description and prediction of the time course 
of drug effects resulting from different drug administra- 
tion regimens. 

Drug concentration-response curve: PD model 

In general, the magnitude of a pharmacological effect 
increases monotonically with increased dose, eventually 
reaching a plateau level where further increase in dose has 
little additional effect [6] . The classic and most commonly 
used concentration-response model is the Hill equation 
[26], also known as the sigmoidal £ m ax model [27] or 
logistic model [28]. The relationship between the concen- 
tration of the drug and its effect is most often nonlinear. 
In this study, we use hybrid systems to approximate PD 
curves. A common method is to replace certain slowly 
changing variables by their piecewise linear approxima- 
tions (see Figure 7). For example, the PD model used in 
our study can approximate the popular sigmoidal £ max 
model (see Figure 8). The £ max model has the general form 
E = V rm\ r m > where £ max is the maximum effect, C is the 
concentration, EC^o is the concentration necessary to pro- 
duce 50% of £ max , and m represents a sigmoidity factor or 
steepness of the curve. 

We assume a threshold of concentration below which 
the drug candidate is ineffective, the minimum effective 
dose (MinED), and another threshold value, called max- 
imum effective dose (MaxED), above which there is no 
clinically significant increase in pharmacological effect in 
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Figure 7 The PD model: concentration-response curve used in 
this study. 



this study. As an example, we use a linear curve to approx- 
imate the concentration-response curve between MinED 
and MaxED. It is assumed that the drug-effect coefficient 
(the drug target is x\) is related to the concentration u 
through a sigmoid function and can be approximated by 
the curve shown in Figure 7. The corresponding relation- 
ship can be expressed as 



q«(u- Of) 
cfM ~ Of) 



(6) 



where q± is the ratio between the drug-effect factor and 
the effective drug concentration u — Of in the linear range. 
This reflects the fact that the drug only starts to take effect 
when its concentration level is above a lower threshold 
Of and its effect saturates when its concentration level 
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Figure 8 Sigmoidal E max model {m = 4), and approximation by 
our PD model. 



exceeds an upper threshold 0™. Note that the sigmoidal 
£ max model can be well approximated by the proposed 
PD model. By taking the derivative of E with respect 
to C and evaluating it at EC50, we obtain the slope as 
q" = ™£C5o ' ^ ne u PP er an d lower bounds should satisfy 
qU(QU _ qU.^ _ £ max A n example of the sigmoidal £ max 
model when m = 4 and our proposed PD model are plot- 
ted together in Figure 8, where the proposed model closely 
resembles the sigmoidal £ max model. 

Periodic drug intake: PK model 

Drug concentration at the effect-site is critical for its phar- 
macological effect. Currently, plasma drug concentrations 
are markers that serve as surrogates for drug concentra- 
tion at the effect-site for beneficial and adverse effects; 
however, markers not grounded on a sound theoretical 
foundation and therapeutic mechanism-based interven- 
tion can limit the usefulness of PK/PD modeling to drug 
development. For example, it has been demonstrated that 
the intracellular PK of a drug is quite different from 
plasma drug concentration [29,30]. As observed in the 
study by Kuh et al. [29], the intracellular concentra- 
tion of a drug will exponentially increase as the drug 
is absorbed after each drug intake. The drug concentra- 
tion may change very slowly (in our model, we approx- 
imate that as a flat curve) when the intracellular and 
extracellular drug concentration approach equilibrium. In 
time, drug concentration will exponentially decrease as 
the rate at which it is eliminated is more than the rate 
at which it enters the effect-site and, as a result, effects 
diminish. 

Based on the study by Kuh et al. [29], a general model for 
drug concentration-time profile is given in Figure 9. Drug 
concentration is plotted on a logarithmic scale against 
time following each periodic drug intake. k a denotes 
the exponential increase quotient; is the exponen- 
tial decrease quotient; r is the interval between each 
drug intake; and pi, P2, and p$ denote the time stayed 
in the increase, equilibrium, and decrease stage, respec- 
tively. Different drugs work in different ways and the 
proposed model is general enough to cover various cases. 
Drug concentration may increase very quickly and, as a 
result, the increase stage may be neglected, or the equi- 
librium stage may be very short and can be ignored 
for simplicity. By adjusting the parameters in the pro- 
posed model, specific drug characteristics can be rep- 
resented. In the case when the proposed model cannot 
approximate a drugs PK profile, extensive simulations 
can be performed based on the drugs actual PK profile. 
In this article, we consider a periodic drug intake sce- 
nario. Specifically, we are interested in investigating and 
comparing the following two potential scenarios: large 
dose with a longer interval versus small dose with a 
shorter interval. 
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T T 
Figure 9 A general drug concentration-time profile. 

Mathematical analysis of drug effect 

In this section, we study the time course of drug effect 
for different dosage and schedule arrangements where the 
drug is designed to repress a "target gene" The case with a 
special PK profile (drug concentration only has exponen- 
tial decay) was analytically studied in our previous work 
[22]. In this study, we extend the analysis considering a 
general PK profile given in Figure 9 and PD model given in 
Figure 8. Closed-form analytical solution is provided and 
simulations are performed to validate theoretical analysis. 
In later sections, we show that the same methodology can 
be applied to interactive genes, where not only will the 
drug affect the gene expression level, but the target gene is 
also coupled with other genes. 

It is assumed that the disease gene has lost part of its 
self-regulation capacity and the dynamical equation of the 
expression level x\ is given by 



*i = Pi ~ Yixi- 



(7) 



There is a steady state x\ = pi/yi in such a system, 
however, if the synthesis rate is much bigger than the self- 
degradation rate, Pi ^> yi, then the gene expression level 
will be too high. A drug is used as a control input to 
repress the target gene expression level. The correspond- 
ing dynamical equation after drug intake is changed to 



where for kxi < t < (k + l)r; and i = 1, 2, . . . denoting 
the index of different dosing regimens, q^ = q™, 6± l = Of, 
Q^i = qu, for any i, since we assume that the same drug 
is taken in different dosage and schedule settings, f,- is the 
highest concentration level reached after taking the drug. 

State-space analysis 

The state-space and a sample trajectory schematic of the 
state (target gene expression and drug concentration level) 
under periodic drug intake are shown in Figure 10. As 
is common in hybrid systems, there are both continuous 
quantitative changes and discrete transitions in our pro- 
posed model. The entire state space may be divided into 
different domains according to the value of discrete state. 
As shown in Figure 10, there are five domains in the state 
space, with D\, D^Ds not being transient domains. The 
figure shows the case when the drug is effective and the 
drug dosage is large enough that & is higher than the 
upper threshold 0". The sample trajectory of the state cor- 
responds to two periods of drug intake (numbers 1-6 cor- 
responding to the junctions of the drug concentration and 
the boundaries of the domains, also marked in Figure 9). 
Another possible scenario is that the drug dosage is not 
large enough that & is between the upper threshold 0" and 
the lower threshold Of. The third scenario is the case that 
$i < Of and the drug is not effective. 

When the state transits in each period under peri- 
odic drug intake, it may pass through different domains 
(depending on the changes of drug concentration along 
time). During the transit time through domains D5 and 
D3, the gene expression level is pushed lower (to the left), 
while the driving strength will depend on the drugs PD 
characteristic. During the transit time through D\, the 
expression level will rise (to the right), since the drug con- 
centration is lower than Of. For the drug to be effective, 
the reduction of the expression level in D5 and D3 has 
to be larger than the increase of the expression level in 
D\. In summary, we should have x\((k + l)r) < xi(kr), 



xi= Pi - yixi - YiXi, (8) 

where y" is the drug-effect factor defined in the previous 
section. After incorporating a drugs PK/PD (Figures 7 and 
9) into our proposed hybrid system model Equation (8), 
considering the scenario that the patient is taking the drug 
periodically, the resulting model is given by 

xi = fix - Y1X1 - q\ l (ui - Of)s+ '(u,i,0f)s~ '(UbO^)xi 
-q^iOf -O^s+iuuOl^Xi, 
Ui (t)= (e ka{t - kTi) - l)s-(f 9 kTi+pi) 

+ SiS~(t,kti +pi +p 2 )s + (t,kTi +p{) 
+ ^e- X ^-kri-Pi-P2) s + {tikTi + pi + p2 ) s -( t , (k + I)!,), 

(9) 




x((k + 2)r) .v((* + l)r) x(kr) x 



Figure 10 The state trajectory schematic (target gene 
expression versus drug concentration) under PK profile 
(Figure 9) assuming dose > Q\. 
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so that after each treatment the expression level %\ will 
decrease. 

State trajectory analysis 

We analyze the drug effect considering the scenario shown 
in Figure 10, where & > 0". The same methodology can 
be applied to a simpler scenario where Of < & < 6". 
We divide the state trajectory in a period kxi < t < 
(k + l)tj into stages a, b, c, d, e, f, and g as marked in 
Figure 10 and examine the drug effect stage-by-stage. The 
time notations used in the derivation are given by 

• t\ : the traveling time from the initial state to the 
boundary between D$ and D\. 

t 2 : from initial state to boundary between D$ and D%. 
£3: from initial state to the end of stage c, £3 = kti +p\, 
£4: time at which the drug concentration starts to 
decrease, £4 = kxi + p\ + p 2 . 
£5: from the initial state to the end of stage e. 
t^: from the initial state to the end of stage f. 



For kii < t < (k + 1)t/, where i is the index for dif- 
ferent dosing regimens, the corresponding equations and 
solutions for each stage are given by: 

• Stage {a) - D x {kx t < t < £1): 



x\ = Pi ~ Yi*i> 
Yi 



xi(kxi) - 



Yi 



-y\{t-kxi) 



Ui (t) = e Xa(t - kTi) - 1. 
• Stage (b) - D 3 (h <t<t 2 ): 

x\ =Pi-iYi + - Of)] xi 
*i(fi)A(fi) 



(10) 



xi(t) = 



A{t) 



+ 



A(t) 4 

A{t) = e -[WWT)-n)t-£^ t -W] 9 
Ui (t) = e ka{t - kTi) - 1. 



(11) 



• Stage (c) -D 5 (t 2 <t<t 3 = kti + pi): 
*i=A-[n+^F-^f)]*i = 
xi(t) = 



Pi 



yi+q u {0 u_ 0 u_ 0 ^ 



+ 



Xi(t 2 ) 



Pi 



Ui(t) 



Yi+q^f) 

x e -[n+fiV?-8T)Kt-t 2 ) 9 

e X a {t-kxi) _ 1 



(12) 



• Stage (d) -D 5 {t 3 <t<t 4 = kx t + pi + p 2 )- 
Xi =fr-[ Yi + q?(6f - Of)] Xl =► 

X\(t) = tt 



+ 



x e -[Yi+q"(of-or)Kt-t3) > 

Ui (t) = w™ x = f, 
• Stage (e) - D 5 (t 4 < t < t 5 ): 

Xi =fr-[ Yi + q u M ~ e f)l *l => 
ft 



(13) 



*i(f) = 



Yi + q\(0"-of) 



+ 



Xi(t 4 ) 



Pi 



Yi + qlW-ov. 

x e -[n+9?(ef-9 1 -)](t-t4) j 

= r/e" Ad(t_i4) - 

• Stage (/•) -£> 3 (*5 < * < fe): 

*i =ft-[yi +q t [(u- 0f)]#i =>■ 



(14) 



xi(f) = 



«i(fe)e La ^ 1 



+(9?ef-n)te] 



Ac 



Ui (t) = 0fe- Xd(t - t5 \ 
• Stage (gf) -Di (£ 6 < £ < (k + l)r,): 



(15) 



xi = Pi- yixi 
Yi 

Ui {t) = 6"e- kd(t - t6 \ 



xi(t) 



xi(t 6 ) 



Yi 



-Yl(t-t 6 ) 



(16) 



We can deduce the necessary and sufficient condition 
for the effectiveness of the drug by expressing the inequal- 
ity xi((k + l)r) < xi(kr) in terms of dosing period r 
and unit dose, assuming the dosage is proportional to the 
maximum drug concentration reached after taking the 
drug. When the initial conditions are x\ = #i(At/), the 
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equations governing the state trajectory from time kit to 
time (k + 1)t; are given by 



xi(ti) = h 



Xi(t 2 ) = 



Yi 

xi(ti)A(ti) 



xiikti) - 



Pi 
Yi 



-Yiih-kxi) 



(17) 



A(t 2 ) A(t 2 ) 



A(t 2 ) = e 
xi(t 5 ) 



x [ 2 Ac-K^^-^-S^^]^, (18) 
Jti 

[(^(l+df)- n ) fe -g^(*2-«)] 



+ 



Pi 



(19) 



p -[n+^(^-^f)](fe-f2) 



* l( ; 5) r [ ^ +( ^-^ ] 



q ±e»e- k d«6-t5)+ {q uef- yi )t e ] 



(20) 



((* + i) T/ ) = ^ + L (*>) - ^1 c -»«*+i)«-«, (21) 

n L n J 



*i=*ti+T- ln(l+0f), 



= Art; + — ln(l + 6(f), 



Pl = —ln(l + fc), 

'-a 



1 



= ts - — In , 



(22) 
(23) 
(24) 

(25) 

(26) 



For the drug to be effective, we need the disease gene 
expression level to decrease following each period of drug 
intake. Hence, we can express x\((k + l)r) < xi(kr) 
in terms of dosage and frequency schedule and derive 
the region where the drug is effective using the above 
listed equations. 



Results and analysis 

Based on the theoretical analysis in previous two sections, 
we demonstrate that the drug efficacy depends on total 
drug intake, different dosages, and frequencies. The den- 
sity of drug intake is denned as a = ^ = It is 
proportional to the total drug intake, and hence, related 
to the drug toxicity level in practice. First, we demonstrate 
the effect of total drug intake (equivalently, a) towards 
drug efficacy. For each fixed total drug intake, we plot the 
target gene expression reduction based on Equations 17 to 
26 for different dosing period r as a curve in Figure 11. It 
is observed that the curve is "U" shaped and there exist 
"sweet spot" for certain dosages and schedules given a 
fixed a. If we define drug efficacy region (DER) as the drug 
drives down the target gene expression by more than a 
desired percentage (say 60% in this case), it is demon- 
strated that the DER is related to the total drug intake, 
dosing period r and dosage. DER is marked by the shaded 
area in Figure 11 for the case that a < 0.5. It is also 
observed that when a gets bigger, which indicating more 
toxicity, DER is getting bigger accordingly. We would like 
to emphasize that toxicity is one of the primary causes 
for drug attrition and long development cycle times [31]. 
If a drugs toxicity profile is available, for example, the 
maximum dosage and maximum exposure (a), we can 
find a good compromise between toxicity and drug effi- 
cacy based on such study, and determine the "sweet spot" 
(a good dosage and schedule balance) given the obtained 
a, and hence provide valuable suggestions of the dosing 
regimens to the clinical study. 

Second, we test the analytical results via numerical sim- 
ulation using MATLAB/SIMULINK. Given a fixed total 
drug intake, or equivalently, a fixed density of drug intake 
(a = 0.4), three typical scenarios are studied by simula- 
tion: small frequent drug intake with r = 7 and dose = 
2.8; big infrequent drug intake with r = 22 and dose = 
8.8; and intermediate dosage and frequency with r = 12 
and dose = 4.8. The results are shown in Figure 12a- 
f, with the first row corresponding to the state responses 
and the second row corresponding to the state space tra- 
jectory. Although the three cases have the same total drug 
intake and initial condition (initial gene expression level 
x(0) = 20), the drug efficacy is different. In the small 
frequent intake case, the dosage is small and the drug con- 
centration is mostly changing between domains D% and 
D\. Figure 12d shows that the state-space trajectory set- 
tles in a small limit cycle and disease gene expression level 
settles at 11.5 at the end of each period of treatment. On 
the other hand, the big infrequent drug intake case results 
in a big limit cycle as in Figure 12f. Although the dosage 
is high, the long period between dosages means that the 
period stayed in D\ is getting longer (where drug con- 
centration is below Of, hence not effective), and disease 
gene expression level settles at 12.1 at the end of each 
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Figure 1 1 The percentage change of the disease gene expression versus the period of drug intake r for a = 0.4, 0.5, 0.6, 0.7, 0.8, 
respectively. The change must be significant (say at least 60%) and a must be smaller than the acceptable toxic level. Parameters used: q\ 



: 0.1, 



% = 1, yi =0.04, e\ 



10, Of 



2, x(0) = 20; p y = 1; p 2 = 5; X d = 0.3. 



period of treatment. As a comparison, the drug effect for 
the case with intermediate dosage and frequency shown 
in Figure 12b,e is superior to the other two cases. Dis- 
ease gene expression level settles at 8.8 at the end of each 
treatment period. If we check the curve in Figure 11 with 
a = 0.4, the intermediate dosage case with r = 12 
is located near the bottom of the "U" shape. Lastly, we 
observe that all state-space trajectories follow the state 
trajectory schematic in Figure 10, as predicted by the 
analytical results. 

Analysis of 2-gene networks 

We extend the theoretical analysis to 2-gene networks and 
show that the same framework applies to the modeling 
and analysis of drug effect in more complex gene net- 
works. We assume that x\ is the target gene, there exists 
a positive feedback loop between x\ and another gene 
X2, and that a drug targets x\ by reducing its expression 
level and providing a negative feedback term —y™x\. The 
resulting 2-gene network is shown in Figure 13 and is 
given by 



X2 = P2 + m*\ 



yixi 

Y2*2 



YiXi 



(27) 



where /3i, /3 2 are synthesis factors, y\ and 72 are degra- 
dation factors, and is a drug-effect factor. r)\ > 0 
and rj2 > 0 are the parameters of the positive feedback 
loop between the two genes. The 2-gene network under 



drug perturbation model, Equation (27), can be rewritten 
as a second-order ODE: 



^i+(ki+Ki M +K2)^i+((Ki+Ki M )K2-^i^2)^i 



The solution of this equation is given by 
Xl{t) ~ \k 1 ^ t + k 2 te Klt ki=k 2 ' 



(28) 



(29) 



where A-i and X 2 are the two eigenvalues of Equation 28 
and k\ and k 2 are parameters depending on the initial con- 
ditions. Letting a = yi + Y\ + Y2> b = (yi + Y\)Yi r l\ r l2 
and d = P1Y2 + forjl* the two eigenvalues are given by 
A. u = - a± Jf-* b . It is easy to verify that a 2 - 4b = 
(Yi + Y\ — Y2) 2 + rj2 > 0. Since a 2 — 4b > 0, we 
conclude that k\ ^ X 2 and both eigenvalues are real. 
Furthermore, one of the eigenvalues, X 2l is always nega- 

-a— Vcfi 



tive since k 2 = ~ a -^~ 4b 
determined by the sign of b: 

X 1 < 0, if b > 0 
= 0, if b = 0 
Ai > 0, if b < 0 

In other words, 



Ai < 0, x 2 < 0 if (71 + Y\)Y2 > mm 
Ai = 0, a 2 < 0 if (71 + Y\)Y2 = mm 
Ai > 0, a 2 < 0 if (71 + Y\)Y2 < mm ' 



0. The sign of Ai will be 



(30) 



(31) 
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Figure 12 Drug response at three different schedules but same drug intake: (a-c) the state response at t = 7, 12, 22, respectively; (d-f) 
the state space trajectory at r = 7, 12, 22, respectively. Other parameters are q u } = 0.1 , fa = 1 , y } = 0.04, Of = 1 0, ef = 2, x(0) = 20; 
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7, 

Figure 13 A 2-gene network with positive feedback loop and 
drug input. 



The above equation has an important biological inter- 
pretation: when the degradation of x\ due to the strength 
of the drug is faster than the increase of x\ due to the 
positive feedback loop, both eigenvalues are negative, the 
system is stable and x\ will experience exponential decay; 
on the other hand, if the effect of the positive feedback 
loop is dominant, then one of the eigenvalues will be 
positive and x\ will increase exponentially. 

Given initial condition x\(to) and x\(to), then for the 
case A.i ^ A.2, we have 

ki = * \° (dX 2 /b + xi (t 0 ) - X 2 xi (t 0 )) (32) 

Al — A2 

k 2 = —(dk^b + Xiito) - kixxfo)) (33) 

A 2 — Ai 

Now with the baseline analysis of the second-order sys- 
tem, we provide detailed state trajectory analysis by taking 
into account the practical form of PK/PD (y") when the 
drug is taken periodically. 

State trajectory analysis 

We analyze the drug-effect following the same framework 
given in the subsection "State trajectory analysis" under 
the main section "Mathematical analysis of drug effect". 
For kii < t < (k + 1)t/, i = 1,2, . . ., the correspond- 
ing equations and solutions for each stages are given as 
follows: 

• Stage (a) - Di (kti <t<t\): 

X\ = + T\\X2 - Y\X\ 

X2 = P2 + mx\ - Y2X2 

Ui (t) = e ka{t - kXi) - 1. (34) 

The solution of x\ (t) is given by Equation (29), with k\ 
and I<2 given by Equations (32) and (33), and to = kx{. 



• Stage (b)-D 3 (ti < t < t 2 ): 

x\ = P\ + mx2 - Yixi - Y\X\ 

*2 = P2 + mx\ - Y2X2 

x\ + ax\ + bx\ — d 

Ui (t) = e Xa{t - kxi) - 1 (35) 

where a, b, d are denned as before. When 
incorporating the practical form of Y\ = q^(u — Of) 
and u = ^{t-kxi) _ ^ fa e ^ove second-order ODE 
has no closed-form solution. In this case, the solution 
can be obtained numerically. 

• Stage (c) - D 5 (t 2 < t < t 3 = kxt + pi): The set of 
equations are the same as in Stage (b) except that 

= q"(0" — Of). Since 7" does not depend on 
u — e Xa ^~' <Ti ^ — 1 explicitly, x\ has a closed-form 
solution given by Equation (29). 

• Stage (d) -D 5 (t 3 <t<t 4 = kti +pi+ p 2 ): The 
solution of x\ is the same as that in Stage (c) except 
the start and end times, and the equation of u, which 
is Ui(t) = u™ x = & in this stage. 

• Stage (e) - D5 (£4 < t < £5): The solution of x\ is the 
same as in Stage (c) except the start and end times, 
and the equation of u, which now is 

Ui (t) = $ ie - x d(t-t4) 4 

• Stage if) -D3 (£5 < t < t^): The solution of x\ is the 
same as in Stage (b) except the start and end times, 
and the equation of u, which now is 

Ui (t) = 0fe- x ^- t5 \ 

• Stage (g) -Di (t 6 <t < (k + l)r/): The solution of x\ 
is the same as in Stage (a) except the start and end 
times, and the equation of u, which now is 

Ui (t) = efe-^-te)^ 

We can deduce the necessary and sufficient condition 
for the effectiveness of the drug by expressing the inequal- 
ity x\((k + l)r) < x\(kx) in terms of dosing period r 
and unit dose. In the 2-gene case, no explicit closed-form 
expression can be deduced for the solutions in stages (b) 
and (f ) and numerical methods have to be applied. How- 
ever, through such analysis, it is observed that the same 
methodology for analyzing drug effect can be extended 
to GRNs with multiple interactive genes, although the 
mathematics involved will become more complicated and 
sometimes numerical methods must be applied when 
there is no closed-form solution. 

Simulation results and analysis 

When drug input is not present, the disease gene expres- 
sion will grow unbounded owing to the positive feedback 
loop between the two genes. Here, we study response of 
the disease gene expression to drug input and compare 
two different schedules for r = 20 and r = 30, keeping 
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Figure 14 Drug response follows two different schedules but same drug intake (a = 0.8): (a-c) the state response, state space trajectory, 
and 3D state-space trajectory at r = 20, respectively; (d-f) the state response, state space trajectory, and 3D state space trajectory at 
t = 30, respectively. Other parameters are q\ = 2,0] = 0.1, ft = 0.1,771 = 1,772 = lyi = 0.2, y 2 = 1,^ = 8,0^ = 3,^(0) = 1,x 2 (0) = 1, 
pi = 1,p 2 = 15,A d = 0.3. 
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a — 0.8 fixed. The response and state trajectories in 2D 
and 3D are given in Figure 14a-f, with the first and second 
rows corresponding to r =20 and r = 30, respectively. 
We observe that both cases have periodic responses, but 
the disease gene expression is much better controlled 
when r = 20. This is because the drug concentration is 
high enough in both cases compared to the threshold 
while the decay of the drug concentration is shorter in the 
case when x = 20. In Figure 14c,f, the 3D state (disease 
gene expression) trajectories show that the trajectory set- 
tles to an inner circle when r = 20, whereas the trajectory 
settles to an outer circle when r = 30. Similar obser- 
vations apply to Figure 14b,e. Note the scale of #-axis of 
Figure 14e is 20 times bigger than that of Figure 14b. 

Extension and discussion 

In previous sections, we have considered the drug-effect 
on one-gene and a 2-gene case. In this section, we will 
consider the drug-effect on a target gene in a more sophis- 
ticated GRN context. 

3-gene network with multiple feedback loops 

Suppose a 3-gene network is given by 

xi = P\s~(x2,0l) - y"x\ + 7)1X3 (36) 
x 2 = P2S~(x lf 0f) -y 2 *2 (37) 

*3 = fas~(x2,0l) - Y?>s~(xi,0l)xz, (38) 

where rj\ is a perturbation parameter (from X3 to x\), 
Pi, @2> k are activation factors, 3/2, 73 are degradation 
factors, and Of, 0\, Of, B\ are threshold values, y" is 
a drug-effect factor. We assume periodic drug intake and 
drug concentration level follows exponential decay dur- 
ing each period, i.e., Ui(t) — £ i e~ kd ( t ~ kTi \ where kxi < 
t < (k + l)Tj. A graphical model of the 3-gene network is 
given in Figure 15. There exist two positive feedback loops 
between x\ and #3. 

When the target gene is in GRN context, not only its 
expression level is related to drug perturbation, but also 
depends on network contexts. Several interesting phe- 
nomena are observed through our simulations study: 

1. Drug response is related to disease stage. Simulations 
are performed with different initial target gene 
expression level (#i(0)). Figure 16a-c shows the 
system responses with x\ (0) = 20, which is not too 
high (corresponding to early disease state). As shown 
in Figure 16a, x\ expression level reduces to the 
range [ 7.7, 8.4] under periodic drug intake, while X2 
and X3, the two other interactive genes settle at 1.0 
and 4.0, respectively. The system reaches a new 
steady state (a semi-stable limit cycle, to be exact), 
with 4 = = h+mMn, ^ = & and 

Y\ Y\ Y 2 




Figure 1 5 The 3-gene GRN model. The solid line is the real 
interaction between genes. The dashed line is derived to show the 
positive feedback loop for certain conditions. 



= £|, where x\ is well controlled. The trajectories 
of x\ vs. u and x\ versus X3 are given in Figure 16b,c, 
respectively. The semi-stable limit cycle is shown in 
Figure 16b. 

System responses with x\ (0) = 40 (corresponding to 
late disease state) are shown in Figure 16d-f for 
comparison. Although the other parameter settings 
are exactly the same, the drug will not repress the 
disease gene x\ (Figure 16d) owing to the interaction 
between the disease gene x\ and gene X3. When 
*i(0) = 20 < Of = 21, Equation (36) becomes 
X3 = fi3S~(x2, #f ) — X3#3> and thus X3 is negative 
regulated by x\ and converge to #3 = ^ . However, 
when initial condition x\ (0) = 40 > Of = 21, 
Equation (36) becomes X3 = P3S~(x2, #f ), and thus 
X3 is positively regulated by X2 and its expression 
level will keep increasing. As a result, x\ will keep 
increasing as well, and a positive feedback loop is 
formed between x\ and X3. This is confirmed by the 
trajectories of x\ versus u and x\ versus X3 given in 
Figure 16e,f, respectively. 
2. Under certain conditions, single drug perturbation 
may not be enough. A drug is usually designed to a 
specific target. In this example, the drug tries to 
provide negative feedback to the regulation of x\ 
(tries to repress x\)\ however, since the target gene is 
interactive (or, in a more general setting, pathways 
have crosstalk), only repressing the target gene (or 
blocking the signal of one pathway) may not prevent 
the target gene from expressing itself through 
interactions with other genes (or through 
inter-connected pathways). In our case, x\ is 
interactive with X3. To continue with previous 
simulation (results shown in Figure 16d-f, we try to 
increase the drug dosage tenfold from u(kr) = 24 to 
u(kr) = 240 with the same dosing period x = 8 
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Figure 16 Simulation results for the 3-gene network under drug perturbation with different initial conditions: (a-c) with initial condition 

(#i(0) = 20), (d-f) with initial condition (#i(0) = 40). Other parameter settings are: x 2 (0) = 0.7, x 3 (0) = 0.5, r = 8, u(kr) = 24, q\ = 0.21, 
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trying to bring down the expression level of %\. 
However, from system responses shown in 
Figure 17a-c, it is observed that the drug is not 
effective although the dosage is increased tenfold. 
One step further, not only we increase dosage to 
u(kz) = 240, but also to increase the dosing 
frequency (dosing period is decreased from r = 8 to 
t = 2), systems responses are shown in 
Figure 18a-d, where Figure 18c shows the left part of 
the trajectory shown in Figure 18b. It can be 



observed that although the drug perturbation is very 
strong, and the drug concentration is always staying 
in domain D5, drug is still not effective. 

From the nonlinear dynamical system perspective, the 
equation X s , = ^ 1+ ^* 3 = ^ 1+?7l f 3 ^ K3 represents a semi- 
stable limit cycle. If the initial condition is from the inside 
of the limit cycle, then the system will converge to the limit 
cycle; however, if the initial condition is from the outside 
of the limit cycle, then the system will diverge from the 
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limit cycle. Such simulation results demonstrate the het- 
erogeneity of the drugs responses due to the nonlineari- 
ties in complex systems, where multiple inputs affect each 
output and the underpinning structure may include par- 
allel, redundant, and feedback loop processes, it is likely 
that some cases will not respond to a single drug pertur- 
bation no matter how strong it is. As a result, innovative 
perturbation methods, such as finding a better target or 
combinatorial therapy, are necessary. 

Simulation of effects of different drugs and a drug 
combination on NF — kB pathway 

In this article, the models and examples are selected 
such that they are mathematically tractable and important 
insights can be obtained, and we can verify the theoretical 
results with simulation results. For large-scale networks 
and multiple drugs/drug targets, the proposed model is 
still applicable; however, analytical results may not be 
attainable even for this simplistic model. In that case, sim- 
ulations can be carried out case-by-case. To illustrate this 
point of view, we carried out a simulation study of the 



NF — kB pathway under two different drugs and each drug 
with different drug targets. 

NF — kB signaling regulates inflammation, cell prolifera- 
tion, and apoptosis by increasing the expression of specific 
cellular genes in response to a variety of extracellular lig- 
ands. How to explore therapeutic strategies to prevent the 
prolonged activation of the NF — kB pathway attracts lots 
of attention [32,33]. An ODE model of the NF - kB path- 
way [34] is adopted and the two drugs under consideration 
are drug X (drug A in [35]) and FDA approved drug pro- 
teasome inhibitor Bortezomib (Velcade) [36]. The detailed 
simulation setup is available in the Appendix and the 
SIMULINK model is given in Figure 19. The specificity of 
some drugs to inhibit several of these components of the 
NF — kB pathway is one of the concerns. For example, the 
proteasome which is responsible for the IfcBa degradation 
has many other important functions. Thus, Bortezomib 
modulates a variety of cellular processes that may con- 
tribute to toxicity if the dosage is too high [36] . Hence, we 
design combination therapy to induce a better effect and 
at the same time to contain toxicity to a certain threshold. 
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Figure 1 8 Simulation results for the 3-gene network under drug perturbation (high dosage and short interval): (a-d) state response, 
trajectory of x\ versus u, detailed initial trajectory of x\ versus u, and trajectory of x\ versus * 3 , respectively. Parameter settings are the 
same with Figure 16 except u(kr) = 240 and r = 2. 



To achieve this, drug X is a protein kinase inhibitor, 
which competitively inhibit IKK with the binding kinet- 
ics the same as that of the natural reaction involving 
NF — kB : IkBcl and IKK [35]. While Bortezomib is 
a proteasome inhibitor that will inhibit IkBcx degrada- 
tion, its effect is adjusted through the parameter setting 
related to individual terms for IkBcx and NF — kB : 
IkBcx molecules rescued from inhibition of IicBa degra- 
dation [35]. We first validate the results in [34,35]. In 
Figure 20a, we show that oscillatory behaviors occur for 
NF — kB pathway with constant stimulus. Under this con- 
stant stimulus, it is observed in Figure 20b,c that only 
very high dose of drug X can effectively block NF — kB 
nuclear translocation. Similar observation is obtained for 
Bortezomib from Figure 20d,e, where low drug dosage 
(65% inhibition) is not effective, while the side effects 



are unacceptable when the drug is effective (95% inhi- 
bition). All the above simulation results are consistent 
with those in [34,35]. In this article, we go a step fur- 
ther and consider the combination of these two drugs. 
It is interesting to see in Figure 20f that some com- 
binations of the drugs with non-overlapping toxicities, 
e.g., combined Bortezomib (65% inhibition) and drug X 
(0.2 /xM), might provide enormous benefit by keeping 
the level of nuclear NF — kB low while having tolerable 
toxicities. 

Conclusions and future work 

This article provides systematic mathematical analysis 
and dynamical modeling of drug effect in the GRN con- 
text, where a drug functions as a control input to reduce 
the elevated target gene expression level. A hybrid systems 



Li et al. EURASIP Journal on Bioinformatics and Systems Biology 201 2, 201 2:1 9 
http://bsb.eurasipjournals.eom/content/2012/1/19 



Page 17 of 20 





. * * A /< >\ ,'• »i 



0.1 
0.09 
0.08 
0.07 
0.06 
0.05 
0.04 
0.03 
0.02 
0.01 
0 



- IkBa 

- IKK 
NFkBn 




l y 1 1 I ' 
/\ / 1 i I I 1 



0 



20 30 40 50 

hours 



0.1 
0.09 
0.08 
0.07 
0.06 
0.05 
0.04 
0.03 
0.02 
0.01 
0 



■ IkBa 

■ IKK 

NFkBn 



0.12 



0.08 




v >j » v » »j w 'j J v ' » v 'J v v i 1 w '» v v 



0.03 



10 



20 30 
hours 



- IkBa 
■ IKK 

NFkBn 



40 



50 



0.1 
0.09 
0.08 
0.07 
0.06 
0.05 
0.04 
0.03 
0.02 
0.01 
0 



0 



IkBa 

IKK 

NFkBn 



20 30 
hours 



Figure 1 9 Simulation results for the NF - kB pathway under various drug perturbations with different drug administration, (a) Effect of 
continuous stimulus, no drug input, (b) Effect of drug X (0.2 /xM). (c) Effect of drug X (1.0/xM). (d) Effect of Bortezomib (65% inhibition), (e) Effect of 
Bortezomib (95% inhibition), (f) Effect of combined Bortezomib (65% inhibition) and drug X (0.2 fjM). The detailed parameters are available in the 
Appendix. 
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model is proposed to study the dynamics of the underlying 
regulatory network under drug perturbation. Drug phar- 
macology information is incorporated into drug thera- 
peutic response modeling to demonstrate the significant 
difference in drug effect for different dosing regimens. 
Considering the complicated nature of gene regulation, 
this study is a small step towards quantitative modeling of 
therapeutic effect. We have kept the examples mathemat- 
ically tractable so that valuable insights and reasonable 
predictions can be obtained from theoretical analysis. 



Compared to our previous work [22], where drug effect 
was only studied for a specific PK profile (drug con- 
centration only has exponential decay stage) when the 
drug is targeted to a single gene, three major exten- 
sions are provided in this article: (i) we provide ana- 
lytical results of drug effect under a very general PK 
profile, where three stages of drug concentration change 
(increase, equilibrium, and decrease) are considered; (ii) 
the proposed methodology is applied to interactive genes 
in a GRN context, with detailed analytical derivations 
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Table 1 Definition of variables 

X] 
X2 
*3 
X 4 

*5 
*6 
*7 
*8 
Xg 



NF-kB 
licBa 

IkBol :NF-kB 
NF - K B n 
kBa n 

kBctn : NF - K B n 
IKK 

IKK : IkBol 

IKK : IkBol : NF — k 



for both one-gene and two-gene cases; and (iii) we per- 
form extensive simulations for a more complicated GRN 
setting and explain several interesting observations due 
to multiple feedback loops and the existence of limit 
cycles. 

It is expected that the theoretical framework proposed 
in this article, when correlated to real biological networks, 
can help improve drug development productivity and 
make drug discovery more systematic. During such pro- 
cess, cross disciplinary effort is indispensable. For exam- 
ple, application of such a framework will require exper- 
iments designed to elucidate model parameters, such as 
protein concentration levels and synthesis and degrada- 
tion speeds. While some parameters may be relatively easy 
to obtain, others may be difficult to get based on current 



techniques and model simplification may be necessary; 
nonetheless, the basic hybrid systems model and the con- 
clusions drawn from it, such as the nature of DERs and 
the role of limit cycles, will remain valid, only their par- 
ticular forms being changed to represent experimental 
instantiation of the model. 

Appendix 

See Tables 1 and 2. 

ODE model of the NF - kB pathway 

The ODE model of the NF — kB pathway is adopted from 
[34,35]. 



dt 



d%2 
dt 



dxs 
dt 

dx/\, 
dt 



= —a^,x\X2 + dq.X's — a^x\x% + + d^x^ 
+ deg^ — k\X\ + /<oi*4 



(39) 



— CL1X2X7 + d\x% — a^x\X2 + d/±x<$ — deg\X2 

- tp\X2 + tp2X5 + S S ynthesis*4(£ ~ t) 



= a$x\X2 — dq.X's — (17X3X7 + d\x^ 

+ fe*6 - deg4X3 

■ k\X\ — a 4X4X5 + ^4^6 — ^01 #4 



(40) 

(41) 
(42) 



Table 2 Parameter values 



Parameter 


Reaction type 


Biochemical reaction 


Value 


Unit 


a 4 


Complex formation 


NF -kB + IkBol -> NF - kB : IkBol 


30 


/xM _1 min -1 


07 


Complex formation 


NF-kB: IkBol + IKK -> NF-kB: IkBol : IKK 


11.1 


/xM _1 min -1 


O] 


Complex formation 


IkBol + IKK -> IkBol : IKK 


1.38 


/xM _1 min -1 


d 4 


Dissociation 


NF -kB + IkBol <- NF - kB : IkBol 


0.03 


min -1 


d, 


Dissociation 


NF-kB: IkBol + IKK <- NF-kB: IkBol : IKK 


0.075 


min -1 


d, 


Dissociation 


IkBol + IKK <- IkBol : IKK 


0.075 


min -1 


deg } 


Degradation 


IkBol -> 0 


0.006 


min -1 


deg 4 


Degradation 


NF-kB: IkBol -> NF - kB 


0.0013 


min -1 


ko] 


Transport 


NF - KBn -> NF-kB 


0.0048 


min -1 


tP2 


Transport 


IkBolh -> IkBol 


0.025 


min -1 


k 2 


Transport 


NF - KBn : kBan -> NF - kB : IkBol 


0.84 


min -1 


ky 


Transport 


NF - kB -> NF - KBn 


5.4 


min -1 


tp1 


Transport 


IkBol -> kBan 


0.05 


min -1 


T 


Synthesis (delay) 


NF - KBn -> NF - KBn + IkBol 


40 


min 


koi 


Inactivation 


IKK -> 0 


0.002 


min -1 


r 4 + d 4 


Catalyzed degradation 


NF-kB: IkBol : IKK -> NF-kB + IKK 


11.1 


min -1 


ri 


Catalyzed degradation 


IkBol : IKK -> IKK 


2.22 


min -1 


■^synthesis 


Synthesis 


NF - KBn -> NF - KBn + IkBol 


0.24 


min -1 
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dxs 

—— = tp\X2 - tp2X$ - a^Xs + (Z4X6 

at 



— a4.X4.Xs — d/^Xfr — /<2#6 



(43) 
(44) 



dX(y 

dt 

dxj 

— = k(t) - k 02 x 7 - aix 2 x 7 + (di + n)# 8 , Ar , 

at (45) 



dxs 
dt 



— a 7 X3X7 + (d\ + Y4)X^ 

a\x 2 X7 — {d\ + r\)x% 



(46) 
(47) 



d%9 , \ 

— — = a-jX^X^ — (d\ + Y4)X^ 

dt 
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